Soft effective interactions between weakly charged polyelectrolyte 

chains 

M. Koniecznyy C. N. Likos, and H. Lowen 
Institut fur Theoretische Physik II, 
Heinrich- Heine- Universitat Diisseldorf, 
Universitatsstrafte 1, D-^0225 Diisseldorf, Germany 
(Dated: February 2, 2008, submitted to The Journal of Chemical Physics) 

Abstract 

We apply extensive Molecular Dynamics simulations and analytical considerations in order to 
study the conformations and the effective interactions between weakly charged, flexible polyelec- 
trolyte chains in salt-free conditions. We focus on charging fractions lying below 20%, for which 
case there is no Manning condensation of counterions and the latter can be thus partitioned in two 
states: those that are trapped within the region of the flexible chain and the ones that are free in 
the solution. We examine the partition of counterions in these two states, the chain sizes and the 
monomer distributions for various chain lengths, finding that the monomer density follows a Gaus- 
sian shape. We calculate the effective interaction between the centers of mass of two interacting 
chains, under the assumption that the chains can be modeled as two overlapping Gaussian charge 
profiles. The analytical calculations are compared with measurements from Molecular Dynamics 
simulations. Good quantitative agreement is found for charging fractions below 10%, where the 
chains assume coil-like configurations, whereas deviations develop for charge fraction of 20%, in 
which case a conformational transition of the chain towards a rodlike configuration starts to take 
place. 

PACS numbers: 82.70.-y, 82.35.Rs, 61.20.-p 
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I. INTRODUCTION 



Polyelectrolytes (PE's) are polymer chains that carry ionizable groups along their 
backbone.i^ Thereby, when dissolved in a polar solvent such as water, these groups disso- 
ciate and a chain that carries charges along its monomer sequence results, leaving behind 
the dissociated counterions in the solution. Upon addition of salt, additional counterions as 
well as coions can be found in the system. One deals therefore with a statistical-mechanical 
problem that entails aspects pertinent to both charge-stabilized colloidal suspensions 4 - and 
polymer solutions. This dual character of polyelectrolytes lies in the heart of the difficulties 
in devising suitable theoretical approaches to tackle questions regarding their conforma- 
tions and structure in nonvanishing concentrations. Indeed, the long-range character of the 
Coulomb interaction renders the typical renormalization-group techniques that have proven 
very fruitful for neutral polymer chains rather inadequate for the problem at hand. More- 
over, the appearance of additional length scales (such as the Bjerrum length, to be defined 
in what follows) has the consequence that several microscopic parameters have an influence 
on the physics of the system, in contrast to the case of neutral polymer chains.— Simulation 
approaches, pioneered by the work of Stevens and Kremer,^ have thus proven very useful in 
shedding light into the question of the typical sizes and conformations of isolated PE chains. 

Considerable work has been devoted to the study of sizes, conformations and charge dis- 
tribution of isolated PE chains under the influence of various physical parameters. These 
include the effects of counterions and salt,— of the temperature,— of annealed versus quenched 
backbone charges 9 as well as of solvent quality . 10 ' 11 ' 12 In the absence of screening, it has 
been found that the typical size R of a PE chain scales with the degree of polymeriza- 
tion N as R ~ A/" (In iV) 1 / 3 at room temperature . 13 ' 14 Another issue concerning isolated PE 
chains that has attracted considerable interest, is that of Manning condensation^ along 
the stiff, rodlike backbone. The fraction of condensed counterions has been examined both 
in theor y 1 ^' 1 ? and in experiments . 1 ^' 19 Counterion condensation and fluctuations have been 
shown to lead to a collapse of the rodlike chains, due to the effective monomer-monomer 
attractions that result from the condensed counterions that position themselves between the 
charged monomers.— ' 21 ' 22 ' 23 ' 24 ' 25 

Relatively less is known about the behavior of polyelectrolyte solutions at finite concen- 
tration. Theoretical approaches on this topic are usually based on integral equation theories 
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that assume rodlike configurations of the chain o 26 ! 27 ! 28 ! 29 ! 30 and make use of the polymer 
reference interaction site model (PRISM) of Schweizer and Curro.— An approach that com- 
bines aspects from field-theoretical methods and liquid-state theory has been developed by 
Yethiraj,— and which allows for a self-consistent calculation of the solvation potential on a 
chain, due to the rest of the solution.— This theory has been applied to calculate the extent 
of the PE chains as well as the correlation functions between the monomers at varying con- 
centration, both ignorin g 32 ! 33 and includin g 34 ! 35 the excluded-volume interactions. Recently, 
phase separation and Coulomb criticality in polyelectrolyte solutions was also studied by 
means of Monte Carlo computer simulations.— 

A very useful theoretical tool that greatly facilitates the analysis of the structure and 
thermodynamics of complex fluids is that of the effective interaction between suitably chosen 
degrees of freedom that characterize the macromolecule as a whole.— For instance, it has been 
recently shown through extensive computer simulations that neutral polymer chains can be 
modeled as ultrasoft colloids if their centers of mass are chosen as the effective coordinates 
that describe the chains in a coarse-grained fashion . 3 ?! 39 ! 4 ^ Effective interactions between 
polyelectrolytes have not been derived up to now, however, with the exception of the case 
of rodlike molecules.— It is the purpose of this paper to fill this gap, by considering weakly 
charged PE chains for which (a) the chains are still flexible and assume a large number 
of conformations and (b) there are no Manning-condensed counterions. We perform this 
coarse-graining procedure by tracing out the monomer and counterion degrees of freedom 
in a salt-free solution and we derive an effective, density-dependent potential of mean force 
between the centers of mass of the PE chains. Since the analytic considerations for the 
derivation of the effective interaction are based on knowledge of the sizes and shapes of 
individual chains, we present first in Sec. |H] simulation and analytical results pertaining to 
single polyelectrolytes and demonstrate the good agreement between the two. The theory 
and the simulation results for the effective potential are presented in Sec. 11111 whereas in 
Sec. HVI we summarize and conclude. Some technical derivations are shown in the Appendix. 
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II. SIZES AND CONFORMATIONS OF ISOLATED PE CHAINS 
A. Simulation model 

We start with the description of our simulation model, valid for both an isolated PE chain 
and two interacting PE chains. We performed monomer-resolved molecular dynamics simu- 
lations (MD). We employed a cubic simulation box with periodic boundary conditions and 
typical edge lengths in the range L = 2R e to 5R e , where R e is the end-to-end distance of the 
chain, depending on the other parameters. To integrate the equations of motion, we used the 
so-called velocity form of Verlet's algorithm . 42 ' 43 ' 44 In order to stabilize the temperature of 
the system, a Langevin thermostat was applied . 45 ' 46 ' 47 This method is based on the introduc- 
tion of convenient random and friction forces taking the fluctuation-dissipation theorem into 
account. The electrostatic Coulomb interaction (see below) of charged particles was treated 
using Lekner's summation method,— where an enhancement of the convergence properties of 
the associated sums is achieved via their adequate re- writing. Thus, a cut-off becomes feasi- 
ble. To limit the memory consumption of the tables containing the corresponding forces, we 
computed the latter by means of a trilinear interpolation," 4 ^ while simultaneously decreasing 
the number of grid points. 

The PE chains were modeled as bead spring chains of N Lennard- Jones (LJ) particles. 
This method was first introduced in investigations of neutral polymer chains and sta ra 6 ' 50 ' 51 
and turned out to be a reasonable approach. To mimic good solvent conditions, we use 
a shifted LJ potential to describe the purely repulsive excluded volume interaction of the 
monomers: 



Here, r is the spatial distance of two interacting particles. <7lj denotes their microscopic 
length scale, and e^j sets the energy scale for the system. We chose the value T = £ljAb 
for the system's temperature, where fee is the Boltzmann constant. 

The bonds of adjacent monomers were depicted via a finite extension nonlinear elastic 
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where &fene is the spring constant. Its value was set to &fene = 7.0£lj- The FENE 
interaction diverges at r = Rq, which determines the maximum relative displacement of two 
neighboring beads. The energy £lj is the same as in Eq. whereas for the length scale 
Rq we have chosen the value Rq = 2.0o"lj- 

Each chain is charged by a fraction a in a periodical manner: every 1/ot bead carries a 
monovalent charge. Due to the requirement of electroneutrality, the same amount of also 
monovalent but oppositely charge ions, namely iV c = aN counterions, are included in the 
simulation box. They are able to freely move in the box, thereby they had to be simulated 
explicitly. 

Finally, the Coulomb interaction Vc ui( r ) between charged particles had to be taken 

into account. Referring to the valencies of monomer ions and counterions as g„ = ±1, 

respectively, the following equation holds: 

Vcoui(^) _ 1 QaQpe 2 
k B T k B T er a/3 

= A B ^, (3) 

where r a/3 = \r a — rg|, with r a and denoting the position vectors of particles a and /3, 
respectively. The Bjerrum length is defined as the length at which the electrostatic energy 
equals the thermal energy: 

Ab = (4) 
ek B T 

where e is the unit charge and e the permittivity of the solvent. For the case of water at 
room temperature one obtains Ab = 7.1 A. In this work, we consider salt-free solutions only. 
Furthermore, the solvent was solely taken into account via its dielectric constant e. The 
Bjerrum length was fixed to Ab = 3.0<tlj. 

The edge length L of the simulation box was set having regard to the simulations' respec- 
tive parameter combinations to suppress surface effects a priori. The center of mass of the 
system was fixed in the geometric middle of the box. The typical time step was At = 0.005r, 



with r = a/ mcTf* j/^lj being the associated time unit and m the monomer mass. The coun- 
terions were taken to have the same mass and size as the charged monomers. 



After a long equilibration time (typically 10 5 time steps), different static quantities were 
calculated during simulation runs lasting about 5 x 10 5 time steps. We carried out simulations 
for charging fractions a = 0.10 to 0.20 and degrees of polymerization iV = 50 to 200. In 
doing so, systematic predictions concerning the a- and iV-dependence of both all relevant 
conformational properties and the effective interaction became possible. 

B. Theory 

In the theoretical investigations of the scaling properties of isolated PE chains, we use 
a mean-field, Flory-type approach, similar to that put forward in Refs. |5 
consider a dilute, salt-free solution of density p c hain = N^^/V, where iV c hain denotes the 
number of chains in the macroscopic volume V. For sake of simplicity, we take all ion species 
to be monovalent. Within the model, a single chain of total charge Q is depicted as a spherical 
object of characteristic spatial extent R. The N c counterions form an oppositely charged 
background, whereas we assume the corresponding density profile to be also spherically 
symmetric. The typical length scale i?w — (4vrp chain /3) -1 ^ 3 of the counterion distribution is 
determined by p c h a in alone. Fig. [T] shows a sketch of the situation, with the shaded domains 
visualizing the relevant length scales of our problem. 

In the following steps, we always restrict our problem to the consideration of a single 
PE chain and its associated counterions. This simplification is convenient, since the N c 
counterions of total charge — Q assure electro-neutrality. Particular attention has to be paid 
to the Manning condensation of counterions on the semi-flexible chains . 1 i 2 i 15 i 17 Let b denote 
the average distance of two adjacent monomer ions along the chain backbone. Condensation 
takes place when the dimensionless fraction Ab/& exceeds unity. The inequality X^/b < 1 
holds for all our parameter combinations and according to Manning theory condensation 
effects do not play any role. It should be a reasonable approximation to neglect them 
completely in what follows. This expectation has been confirmed by our simulations. 

Our main goal is to obtain theoretical predictions for the equilibrium expectation value 
of R. It is determined through minimization of a variational free energy, which reads as 

^(R; R w ) = U$ + F d + F int + S®. (5) 

Here, ug* is a mean-field electrostatic contribution, F e \ and F int represent the influences 
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of conformational entropy and self avoidance of the monomers and Sc 1 ^ is the entropic 
contribution of the counterions. All terms will be described in more detail in what follows. 

The electrostatic mean-field energy Ujp is given by a Hartree-type expression of the 
general form 

^ = ^y dVdV £!W), (6) 

with the local charge density q^(t) to be defined below. On purely dimensional grounds, 
we expect a result that writes as 

PV$> = ^ • *« (£) , (7) 

with (3 = denoting the inverse temperature. The shape of the dimensionless function 

depends on the underlying charge distribution £>W(r) alone. 
The elastic contribution F c \ to the free energy is entropic in nature, written as*& 

(3F el = flF(0) + m?t (8) 

with the equilibrium bond length a and the unimportant constant F(0). It stems from 
an ideal, i.e., Gaussian approximation of the conformational entropy of the chain. For 
the additional, non-electrostatic contribution, F int , of the PE chain, arising through self- 
avoidance, we employ the Flory-type expression^"' 



N 2 

(3F int (R) ~ v —. (9) 

Here, v is the so-called excluded volume parameter, whose value has been discussed 
frequently^™ The monomer volume turned out as a good approximation for the case 
of neutral polymer chains. Without Manning condensation, this guess remains valid for PE 
chains. Thus we chose Vq ~ where we identified the typical monomer length with o"lj. 
Finally, the term is an ideal entropic contribution of the form 

(3S^ = J d 3 r p c (r) [In (p c (r)al) - l] + 3N C In ( , (1 ) 

where A is the thermal de-Broglie wavelength of the counterions. Since Sd is independent 
of the variational parameter R, it contributes a trivial constant only and will we omitted in 
all further steps. It was solely included in Eq. (j3J) for the sake of completeness. 

In order to carry out calculations following Eqs. © to Q, assumptions concerning the 
number densities of all ionic particle species are required. On this basis, the total charge 



density g^(r) can be computed. Due to the small charging fractions a considered here 
and according to polymer theory ffi^iffl it is likely to expect a Gaussian shape of the local 
monomer density p m ( r )- For the counterions, an adequate approximation remains at first 
unknown. MD simulations can be helpful to estimate the local densities. Fig. |21 shows cor- 
responding data confirming the above expectations concerning the monomer density profile 
and proving the absence of any condensation effects. The results regarding the counterion 
density are ambiguous. Admittedly, the concrete shape of the latter distribution is of minor 
importance for several reasons: in general, the local arrangement of particles is of secondary 
relevance when dealing with long-range forces. Since we investigate dilute solutions only, 
the counterion density remains low anyway. Moreover, without Manning condensation the 
effective charge of a PE chain equals the bare charge and we expect the monomer ions' con- 
tribution to electrostatic energy to be dominant, while the influence of the counterions 
should be weak. Thus, our main guideline should be to preserve mathematical simplicity. 
In what follows, we will consider two varying sets of density profiles. 



1. Model A: Gaussian monomer density, homogeneous counterion background 



Within the scope of a first approach, we postulate a Gaussian monomer density profile 
according to our MD simulation results 

N 



Pm(r) 



-{v/Rf 



(11) 



and furthermore assume the counterions to constitute a homogeneously charged background 
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where 0(x) is the Heaviside step function. This leads to a total charge density ^(r), 
written as 



g (1) (r) 

Q 
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Therewith and following Eqs. © and (J7J), we compute the dimensionless function 



S 



It reads as 




Summing up all contributions from Eq. © yields the total free energy J 7 ^' of the model 
system. The value R determining the typical width of the spherical monomer cloud depicting 
the PE chain is then found by numerical minimization. It acquires an explicit density 
dependence through Rw, as usual in charged systems. 



2. Model B: Gaussian monomer and counterion densities 



We now introduce a different modeling of the considered system by replacing the homo- 
geneous counterion background of Sec. Ill B II also with a Gaussian density profile, according 
to the equation 

Pc(r) = -£hre-^)\ (15) 



^' 2 R W 



Regarding the monomer distribution, we still follow the assumptions of the previous section 
and define the local density p m (r) using Eq. (JTTJ). For the overall charge density ^(r) we 
obtain analogous to Eq. (JTJJ): 
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In anticipation of our considerations for two interacting PE chains [see Eq. 
IIII Aj . the following equation must hold: 



U§\R; R w ) = ~ Urn Ug>(D; R; R w ), 



r(2), 



(16) 
in Sec. 

(17) 



where U^(D; R; i? w ) denotes the electrostatic mean-field contribution to the free energy of 
two PE chains with center-to-center separation D. Using the results of the Appendix, one 
finally gets 

tf (1) (x) = 



'2tt 



l+x- 



(;)' 



(18) 



Compared to model A, all steps to come are completely analogous. Once more, mini- 
mization results obtained this way exhibit an explicit density dependence via i?w 
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TABLE I: Comparison of the chain radii of isolated PE chains obtained from both theory and 
MD simulation for different monomer numbers N and charging fractions a. The last column 
lists the number N* of free counterions, necessary to calculate the effective chain interaction at 
nonoverlapping distances, see Eqs. (|27jl and (|28j) . 
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C. Comparison to MD result 
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Based on the mean-field models presented above, predictions of the spatial extent of 
isolated PE chains become feasible. In order to verify this results, we have to test them 
against MD data and alternative theoretical approaches. In doing so, we first have to 
address a somewhat technical question. MD simulations yield expectation values for the 
radii of gyration R g and the end-to-end distances R e , but neither of this quantities strictly 
corresponds to the variational parameter R as obtained by minimization of the free energy 
JFW. Consequently, it remains unclear how a comparison can be carried out. Since we 
depicted the PE chains as spherical objects, one should expect R to match half the end-to- 
end distance R e /2. Our results confirm such guess a posteriori, as will be seen shortly. 

Fig. El shows typical MD snapshots, that give rise to some first conclusions. A variation 
of the lateral chain extent becomes manifest, i.e., PE chains are more strongly stretched in 
the mid region than at their ends. This finding is in agreement with the results in Refs. 
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and|lj. In general, the spatial configurations of PE chains are determined by two competing 
contributions to the free energy: the Coulomb repulsion and conformational entropy. For 
monomer ions in the middle of the chain sequence, electrostatic influences dominate and 
the described stretching is energetically favorable. But regarding the outermost parts of 
the chain backbone, the interaction of charged particles is of less importance, while the 
conformational entropy now becomes deciding. This leads to more coil-like arrangements 
of the affected monomers. In particular, short chains (monomer number iV < 50) remain 
coiled even for increasing charging fractions a. 

While simplified theoretical approaches predict a linear relationship R ~ N between 
chain radius and monomer number,— the effect described above requires some modifications 
of this scaling law. Considerations based on the concept of so-called electrostatic blobs yield 
weak logarithmic corrections, resulting in an expression of the general f or m 9 i 13 i 14 

R ~ N(\nN) 1/3 . (19) 

Fig. 0] compares MD data with chain radii obtained by both our mean-field theories and 
the scaling law pursuant to Eq. (|19|) for two different charging fractions. In the latter case, 
the unknown constant of proportionality was conveniently determined by means of a fit to the 
MD results. Thus, theoretical predictions and simulational findings are in good qualitative 
and quantitative agreement. This agreement particularly improves with increasing degree of 
polymerization N, because the basic assumption of continuous density profiles in place of an 
explicit consideration of discrete particles becomes more and more reasonable under these 
circumstances. Moreover, the presented results confirm our identification of the variational 
parameter R with half the end-to-end distance R e /2. Here, we want to point out that the 
radius of gyration is noticeably smaller (see Tab. QJ). 

In addition, Figs. 0] and directly offer a possibility to read off information concerning 
the charge dependence of the spatial conformations of PE chains. According to this, in the 
event of very short chains, the chain radii are approximately unaffected by the charging 
fraction a. Effects due to the Coulomb repulsion of equally charged beads do not clearly 
manifest themselves until the monomer number increases noticeably. Then, with increasing 
a, a transition from coil-like to rod-like configurations takes place. In other words, for fixed 
N, the spatial extent of such chains is monotonously ascending with the charging fraction. 
Deviations from this behavior are likely to expect for highly charged chains only. Under 
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these conditions, Manning condensation becomes more relevant and causes a collapse back 
to coil-like conformations . 1 i 2 i 25 i 61 However, for the range of parameters considered here, we 
do not have to take such processes into account, as our simulations unambiguously prove. 

III. EFFECTIVE INTERACTIONS BETWEEN PE CHAINS 
A. Theory 

To derive a theoretical model for the effective interaction V e s(D) of two PE chains, we 
separately consider the cases of small (D < Do) and large (D > Dq) distances of their cen- 
ters of mass, respectively. The approach we apply is motivated by similar investigations for 
polymer and PE stars , 54 > 62 where an overlap condition for two stars of radius -R s tar determines 
the parameter Dq = 2i? star - But here, the linear chains we deal with are fractal objects of 
a certain flexibility. Thus, the notion of chain overlap in the previous sense becomes mean- 
ingless. The latter fact strongly affects our modeling. According to this, the requirement of 
continuity of the effective potential and the corresponding force when matching the partial 
solutions V^(D), valid for D < D , and Vjf(D), valid for D > D , defines the value of D . 
In particular, D varies depending on the set of parameters considered. 

At first, we assume distances D < Dq. The effective interaction V^(D) results after 
taking a canonical trace over all but the chain center of mass degrees of freedom. Let 
T^ 2 \z) be the Helmholtz free energy of a systems containing two PE chains at center-to- 
center separation z, therewith follows 37 

V cS (D)=^ 2 \D)-^ 2 \oo). (20) 

The term J-'^(oo) is manifestly independent of D and contributes an additive constant 
only. On this account, it does not affect the effective forces, which are obviously obtained 
by deriving the effective potential with respect to D, i.e., using the equation F eS (D) = 
—dV^/dD. Hence, we drop it in all steps to come. 

To calculate the Helmholtz free energy J-'^(D), we start from a mean- field approach 
analogous to Sec. Ill B 21 and model the system of interacting PE chains by superposing two 
charge density profiles Q^(r) as given via Eq. (|TH|) . In other words, we assume them to 
maintain the same profiles as in the case of isolated chains and to completely interpenetrate, 
whereby both charge distributions are shifted with respect to each other by the vector D 
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connecting the chains' centers of mass (see Fig- • Moreover, we expect the parameter R to 
be independent of the center-to-center separation D. Again, the characteristic length scale 
i?w is determined by the density p c hain directly, i.e., by the equation i? w = (47rp chain /3) -1 ^ 3 . 
As can be seen in Fig. the radii of the chains remain essentially unchanged even at full 
overlap. Moreover, in Fig. |S] we show the monomer density profiles, averaged over both 
chains. These are measured during the MD runs as a function of the distance from their 
common center of mass, for various separations D between the centers of mass of the two 
chains. It can bee seen that for distances D as small as IOulj, the individual profiles can 
still be described by Gaussian functions. Consequently, the overall charge density £>^(r; D) 
reads as 



By now, we are interested in the effective interaction as a function of the center-to-center 
separation D exclusively, whereas conformational properties should not be considered. Thus, 
contrary to Sec. Ill Bt we omit any minimization of the free energy. Instead, we use the 
width R of the monomer density profiles as additional fit parameter. As we will see, optimal 
agreement between theoretical predictions and MD data is achieved when R approximately 
matches half the end-to-end distance i? e /2 of isolated PE chains (see Tab. ITT]). 

In what follows, we completely neglect the conformational entropy F e \ and excluded 
volume interaction F int , which are only weakly dependent on the inter-chain distance D. 
Hence, they yield constant contributions only and do not influence the effective potential 
Ves(D) significantly. Our problem immediately reduces to the computation of electrostatic 
term (D) and counterion entropy Sc(D), i.e., we have 



Taking Eqs. (fT7)J) and (j2H) into account, we are able to conveniently rewrite Eq. ©. 
Dropping all irrelevant .D-independent self energy terms, we obtain 



On purely dimensional grounds, the result of the above integration must read as (cf. Sec. 




(21) 




(22) 




(23) 





(24) 
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Once more, the specific shape of the dimensionless function (x) depends on the un- 
derlying charge densities g^ 1 ' (r) and g^ 2 ' (r) alone. A detailed derivation in case of Gaussian 
density profiles for each chain and for the associated counterions (Model B) according to 
Eq. (jl6j) of Sec. Ill B 21 is presented in the Appendix. The corresponding result reads as 



y D 7 D J \2D 2 J \2D 2 J V AD 
where we introduced the abbreviation 

2 r 



h(x) = - dt — ^e"' x . (26) 

7T J Q t 



Since an analytical computation of both h(x) and the counterion entropy S C (D) is not 
feasible, we need to make use of numerical methods, 4 ^ starting from Eqs. (|10|). ()16jl and 
()2ip. By means of this, the Helmholtz free energy Ti{T)\ R; i?w) an d the effective potential 
V^(D) are known in principle. They exhibit a typical explicit density dependence via i?w 
The corresponding forces must inherit such property, what has to be regarded during any 
comparison to MD data. 

Fig. El quantitatively compares the parts Ujp(D) and S$?\D) for an exemplary chosen 
parameter combination. It can be seen that their variations with D, which determines the 
effective force through the D-derivative of these contributions, are comparable, i.e., neither 
of them dominates the physics of our system over the other. This is contrary to comparable 
results for PE stars,— where the electrostatic energy constitutes the major influence for 
small arm numbers /, while it is of minor importance for the opposite case of big /. For the 
latter situation, the PE stars nearly act as if the were neutral, due to the strong adsorption 
of counterions in their interior , 53 i 54 a mechanism absent in the case of PE chains at hand. 

Now, the case of inter-chain distances D > Dq has to be considered. We assume the 
chains to be spheres of average radius R g , whose respective net charges Q* = N*e are 
reduced compared to the bare charges Q = N c e due to counterions located within this 
spheres. Each chain is surrounded by a cloud of the remaining N* free counterions, which 
cause an additional electrostatic screening of the interaction. Again, i? w = (47rp c h a in/3) -1 ^ 3 
denotes the clouds' characteristic spatial extent. 

Adopting a Debye-Hiickel approach 4 for nonoverlapping polyelectrolytes, we postulate a 
Yukawa- type expression for the effective interaction potential there, i.e., 

P -D/X 

PV+{D) = N 2 Xb • (27) 
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Here, the screening length Ad is mainly determined by the radii R g of the spheres and the 
number N* of free counterions. We obtain 37 



F?3 _ p3 

W g 

3N*X B 



1/2 



(28) 



In what follows, we fix both R g and N* according to MD results (see Tab. HJ. Thereby, 
the latter is determined by calculating the time average of the number of counterions within 
an imaginary sphere around a single chain's center of mass whose radius equals the instan- 
taneous value of the corresponding radius of gyration. In doing so, we know the mean net 
charge Q* at the same time. There are no further fit parameters, i.e., the Yukawa tail V^(D) 
of the effective potential is completely specified by Eqs. (|27jl and (|2*8j) . It exhibits explicit 
and implicit density dependences via the quantities i?w and N*, respectively. 

In order to determine the interval boundary Dq, we simply introduce a continuity con- 
straint for the effective force F e g(D) = —dV e s(D)/dD. In mathematical formulation, the 
corresponding condition reads as 

d[V«(D)-y+(Dj\ 



3D 



0. 



(29) 



Do 



Vefl(D) 



The final step in deriving the total potential V g r (D) is to match both parts of the function 
at D — D . This is feasible by constantly shifting V~^(D) only, since V e g(D) must vanish 
for increasing chain separations D —>■ oo. Thus, the following definition is convenient: 

V cS (D) + V+(D )-V cS (D ), D<D 

V£{D), D>D 
Fig. El shows corresponding results for several parameter combinations. The effective 
interaction is purely repulsive in nature, ultra-soft and even bounded for vanishing distances 
D — > 0. There is a strong charge dependence of the potential, i.e., changing the charging 
fraction from a = 0.10 to a = 0.20 roughly yields a factor two. To emphasize the latter 
fact, Fig. HUT c) additionally includes an effective potential for neutral polymer chains, 
which was computed using the equation^ 

D 

_1.13-.R g 

Compared to charged systems, the neutral polymer effective interaction is about two orders 
of magnitude weaker. Within the scope of our considerations, we did not investigate the 
explicit and implicit density dependences of the effective interaction potential in more detail. 



PV$(D) = 1.87 -exp 



(31) 
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B. Comparison to MD results 



Now, our theoretical model has to be tested against MD results. A straightforward 
comparison is not feasible, since (effective) interaction potentials are not accessible directly 
using standard (MD) simulation techniques. By means of the expressions 

F L m) ( r li ' r 2io) = -V rcti() Kff(riio) r 2io) 

= F^(D) (32) 

and 

Fi c )(R 1; R 2 ) = 




(33) 

respectively, we can measure the mean forces acting on given monomers io or the centers 
of mass. Here, let f Q j be the instantaneous force that monomer % of chain a experiences 
and r ai the corresponding spatial position. In addition, the vectors R a denote the chains' 
respective centers of mass. In both equations above, (. . .) D denotes an averaging with fixed 
chain distance D = |r lio — r 2 j | or D = |Rx — R 2 |. For symmetry reasons, the following 
relation must hold: 

While the validity of Eq. is manifest, Eq. (|3*3^1 needs a mathematical proof, described 



in detail in Ref. 



63l In all further steps, we consider the magnitude of the force. It is 



connected to the potential via 



/•;,r(£>)= F$ C \D) 

dV e fi 



(35) 



dD 

According to this, predictions for the effective forces become computable starting from 
theoretical results for the corresponding potential. This allows the desired comparison to 
MD data in principle. 
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Choosing the centers of mass as effective coordinates is arbitrary; one could have picked, 
e.g., the central monomers as representatives of the whole chain. The concrete choice of 
effective coordinates does affect the effective force, as physically expected. However, once 
the separation between the centers of mass or the central monomers exceeds the typical 
chain size, it is reasonable to expect that it should be largely irrelevant which degrees of 
freedom are kept fixed, since the chains appear at these scales as diffuse objects. The latter 
assumption has been confirmed by our MD simulations. Following this argument, we limit 
ourselves to the consideration of the centers of mass as effective coordinates according to 
Eq. flB). 

The mean-field model presented in Sec. 1111 Al predicts ultra-soft and even bounded ef- 
fective potentials V e ^(D) (Fig. ITU]) . Thus, the corresponding forces reach their respective 
maxima at inter-chain distances D/R g > 0, which is in good quantitative agreement with 
MD results. Fig. 1121 shows a concluding comparison of theoretical results and MD data for 
various parameter combinations. For the lower charge fraction, a = 0.10, we find good agree- 
ment between theory and simulation, although we observe qualitative deviations for very 
small inter-chain distances. There, the theoretically predicted forces are in part stronger 
and, in particular, do not drop to zero for D/R g — > 0, as one would expect on symmetry 
grounds. This discrepancy is due to the use of Gaussian density profiles (cf. Sec. lIII A"]) , since 
such approximation becomes unreasonable for vanishing chain separations (see Fig. [SJ). The 
quality of the theory decreases for the larger value of charge fraction, a = 0.20, since the 
chains in this case become more stretched and the underlying theoretical assumption of a 
Gaussian monomer profile starts losing its validity. Here, modeling the chains as rods would 
probably yield better results, although these rods are not completely stiff, i.e., pronounced 
lateral fluctuations are still present. For all plots, values of the fit parameter R (in case 
D < Dq) are explicitly given. They roughly equal the half of the end-to-end distances R e of 
the chains, as measured within simulations (cf. Tab. H|. This again allows an interpretation 
of R as spatial extent of the PE chains and establishes the consistency of the approach for 
the effective interaction with that for the isolated chains. 

The interval boundary D is identifiable by means of a slight cusp in the force vs. distance 
curves. The latter is not physically reasonable, but an artifact arising as consequence of the 
matching of both branches. Without the introduction of at least one additional fit parameter, 
such effect is not avoidable. Our results exhibit a distinct charge dependence of Dq. While 
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TABLE II: Comparison of the fit parameter R to the chain radii of isolated PE chains as obtained 
by MD simulation for different monomer numbers N and charging fractions a. 



N 


a 


(R/a L] ) a 


(R s /a L i) b 


(R e /2a L] ) b 


50 


0.10 


6.1 


5.7 


7.5 


50 


0.20 


9.0 


6.6 


9.2 


100 


0.10 


12.7 


10.4 


13.6 


100 


0.20 


20.5 


12.4 


17.4 


150 


0.10 


19.4 


13.9 


21.0 


150 


0.20 


31.4 


21.0 


32.8 


200 


0.10 


24.7 


23.7 


33.0 


200 


0.20 


44.0 


28.9 


42.8 



a Theory. 

''MD simulation, iV c h a in = 1- 

Do/Rg « 2 holds for a = 0.10 and independent of the monomer number N, the position of 
the matching moves towards bigger distances D /R g ^4 for a = 0.20 and arbitrary N. 

IV. SUMMARY AND CONCLUSIONS 

We have presented a theoretical approach in describing the conformations and sizes of 
polyelectrolyte chains. The theory is based on the assumption of Gaussian monomer- and 
counterion-profiles around the chain's center of mass and employs a variational free energy 
that includes electrostatic, excluded-volume and entropic contributions. By direct compar- 
ison with simulation results, we have shown that the theory is capable of predicting the 
typical size of PE chains that are weakly charged and its dependence on the charge fraction 
and the degree of polymerization. Thereafter, the theory has been extended to describe 
two interacting polyelectrolytes and the effective interaction potential between their centers 
of mass has been derived and compared successfully with computer simulation results. In 
doing so, we have made a first step in describing polyelectrolytes as soft colloids, in analogy 
with recently developed approaches for neutral polymers .^^2^°. The effective potentials ob- 
tained are ultrasoft and bounded but nevertheless much more repulsive than those obtained 
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for polymer chains. The physical reasons lie both in the electrostatic repulsion between 
the charges carried on the chains and on the entropically caused osmotic pressure of the 
counterions that are trapped within the interior of the chains. 

The effective potential includes an explicit density dependence arising from the redistribu- 
tion of counterions inside and outside the chains upon a change of the overall concentration. 
In principle, this potential can be employed in order to study the structural characteristics of 
concentrated solutions of polyelectrolytes, i.e., the correlations between the centers of mass 
of the chains. Moreover, it can be used for the calculation of thermodynamic properties 
of concentrated solutions, such as free energies and pressure and of phase transitions, such 
as the phase separation investigated recently by computer simulations.— Indeed, whereas 
it has been already shown that the approach of coarse-graining successfully predicts phase 
separation in polymer blends,— the question of whether the same is true for polyelectrolytes 
has not been examined to date. We plan to return to this problem in the future. 
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APPENDIX: PROOF OF EQ. flUD 

In Sees. |n] and IIH[ we dropped technical details of the calculation of the Helmholtz free 
energies. Here, we want to carry out the derivation of the mean-field contribution U^\d) 
for interacting PE chains (cf. Sec. 1111 Aj) explicitly. We start from the equation 



Using f!21l) in (jA.l)) and neglecting all irrelevant self-energy terms, which are independent of 
D, we get 




'( 2 )(r;D)f?( 2 )(r';D) 



(A.l) 




(A.2) 
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Let * denote a convolution of two functions. With the definition v(r) = l/(er), we rewrite 
(jA.2j) obtaining 

U$\d) = j dV q {1 \v' - D) J d 3 r v{\r - v'\)q^\v) 
= J dV £? {1) (r'-D) [f? (1) *v] (r) 

= [f? (1) * [f? (1) *v]] (D). (A.3) 

Now, we pass into Fourier space and thereby introduce the notation / for the associated 
Fourier transform of a function /. First of all, this yields the relation 

[ g W * v] (r) = — ^ [ d 3 k £ {1) (k)?)(k)e ikr . (A.4) 
(2tt) J 

Taking (jA.4|) into account, ()A.3|) reads as 

U h( D ) = [q {1) * * v ]] ( D ) 

d 3 fc ^ (1) (k)[f ? '«Tt;](k)e ik - D 



" (2tt) 3 

With Q^ l \v)/Q = p m (r)/N — p c (r)/N c and regarding the linearity of the Fourier transform, 
Eq. (jA.5|) writes as 

US\D) =-^3 Jd*k p 2 m (k)S(k)e ik ' D + ± Jd*k pt(k)v(V^ 



2 



NN C I ^ P m ( k )Pc(k)^(k)e ikD } . (A.6) 

Up to now, we did not assume any specific properties of the density profiles or the integration 
kernel v. All considerations are valid in a very general fashion. In our special case, the 
mentioned functions are radially symmetric. Due to this fact, we obtain 

U " iD)= ^D\N!J dkksm(kD)p 2 m (k)v(k) + W2 J Q dkksm(kD)pl(k)v(k) 

dk ksm(kD)p m (k)p c (k)v(k) 



NN cJ0 



Q ~ {h(D) + I 2 (D)-2I 3 (D)} (A.7) 



2n 2 D 
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Starting from Eqs. (|16|). (J2~T|) and using the above definition of v, we find 



Pm(fc) _ 
AT 

Pc(fc) 

47T 

5 « - 



-fc 2 B 2 /4 



-fc 2 i? 2 /4 



(A.8) 
(A.9) 
(A.10) 



Based on (|A.7jl and including Eqs. (|A.8|) to (jA.lOjl . we firstly have as an intermediate result 

4tt r °° 



h(D) 



Ait 
2ir 2 



djfc sin(fcD) c _ fc2fl2/2 



t 



■ h 



R 2 



e \2D 2 
In complete analogy, we achieve the relations 



h{D) 



2^ 2 h (K 

1 2D 2 



and moreover 



h(D) 



2rr (R 2 + R 2 N 



e V ^D 2 

For simplicity reasons, we thereby introduced the abbreviated notation 



2 roo 

h(x) = — 

n .'0 



sin(t) _ t2;r 
dt — — e . 
t 



Furthermore, we define based on the previous Eqs. (|A.11|) to (jA.13|) : 



(A.ll) 



(A.12) 



(A.I3) 



(A. 14) 



0(2) / R Rw 

1 £>' D 



2tt 



2 {h(D) + I 2 {D) - 2h{D)} 



h 



R 2 



P 2 
h' W 



2D 2 J \2D 2 
Substituting (jA.15|) in ()A.7|) finally yields the result 



2h 



R 2 + Ran 
AD 2 



R Ryj 



eD 

N c X * . 0(2) 



D 



D' D 



(A.15) 



(A.16) 
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FIGURE CAPTIONS 



FIG. 1: A sketch showing a PE chain of typical spatial extent R and the surrounding 
counterion sphere of radius i?w 

FIG. 2: Local number densities of (a) the monomers and (b) of the counterions for an 
isolated PE chain with N = 200, a = 0.10 and i? w /o"Lj = 124.1, as obtained by MD 
simulations. Here r denotes the distance from the center of mass of the chain. 

FIG. 3: MD snapshots of isolated PE chains with degree of polymerization N = 200 and 
two different charging fractions a = 0.10 (left) and a = 0.20 (right). Bright gray balls are 
neutral monomers, and the dark spheres along the chain indicate the charged monomers. 
The counterions are the small, dark gray spheres around the chain. In addition, the black 
crosses mark the centers of mass of the instantaneous configurations. 

FIG. 4: Spatial extent of isolated PE chains as a function of the degree of polymerization 
N for charging fractions: (a) a = 0.10 and (b) a = 0.20. The lines show theoretical results 
while the symbols denote data obtained by our MD simulations. 

FIG. 5: Chain radii of isolated PE chains as a function of the charging fraction a for 
fixed monomer number iV = 200 and -Rw/°"lj = 124.1, as obtained by our MD simulations. 

FIG. 6: A sketch of two overlapping PE chains, for a more detailed explanation of the 
theoretical modeling, see the text. 

FIG. 7: Simulation results for the radii of interacting PE chains, depending on their center 
of mass separation D. The results shown here refer to parameters N = 100 and a = 0.10. 
The corresponding radii of isolated chains are shown as horizontal lines for comparison. 

FIG. 8: Density profiles (a) for the monomers and (b) for the counterions of interacting 
PE chains with N = 200, a = 0.10, -Rw/clj = 157.6 and varying center of mass separation 
D, as obtained via MD simulation. 

(2) (2) 

FIG. 9: Comparison of the terms (D) and 5c (D) that contribute to the Helmholtz 
free energy as functions of the chain separation D, where we exemplarily chose N = 200, 
a = 0.10 and R w /(?lj = 157.6. 

FIG. 10: Theoretically predicted effective potential V cS (D) as a function of the interchain 
separation D for (a) iV = 100, (b) N = 150, and (c) N = 200. The corresponding charging 
fractions are a = 0.10 and a = 0.20, as indicated in the legends. For all cases, the used 
values of the fit parameter R are specified in the legend boxes. In panel (c) we additionally 
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show as the dash-dotted line the result valid for neutral polymers.— 

FIG. 11: Simulation snapshots of interacting PE chains with degree of polymerization 
iV = 200 and charging fraction a = 0.20. Bright gray balls are neutral monomers, and 
the dark spheres along the chains indicate the charged monomers. The counterions are the 
small, dark gray spheres around the chains. The distances of centers of mass, which are 
marked by the black crosses, are D = 2olj (left) and D = 40<tlj (right), respectively. In 
general, the mean chain directions are not perpendicular to the vector D connecting the 
centers of mass. 

FIG. 12: The effective force F c g(D) on the center of mass of a polyelectrolyte as a 
function of the interchain separation D, and for charging fractions a = 0.10 [(a), (b), (c)] 
and a = 0.20 [(d), (e), (f)], where the respective degrees of polymerization are: N = 100 
[(a), (d)], N = 150 [(b), (e)], and iV = 200 [(c), (f)]. Solid lines are theoretical predictions, 
with the fit parameters R specified in the legend boxes. Data points denote MD results for 
the chains' centers of mass as reference points. 
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FIG. 1: Konieczny, Likos, Lowen 
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FIG. 2: Konieczny, Likos, Lowen 
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FIG. 3: Konieczny, Likos, Lowen 
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FIG. 4: Konieczny, Likos, Lowen 
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FIG. 5: Konieczny, Likos, Lowen 
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FIG. 7: Konieczny, Likos, Lowen 
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FIG. 8: Konieczny, Likos, Lowen 
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FIG. 10: Konieczny, Likos, Lowen 
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FIG. 11: Konieczny, Likos, Lowen 
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FIG. 12: Konieczny, Likos, Lowen 
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